# aem.single.pond -- 
#
# Copyright (c) 2000, V.Kelson and M.Bakker

from types import *
from cmath import log,pi
from math import sqrt

import geometry
XY = geometry.XY

import base

import aem.container
RowList = aem.container.RowList
ColumnList = aem.container.ColumnList
TraceCheck = aem.container.TraceCheck
 
class Pond(geometry.two_d.Circle,base.Element):
	'''Class for circular ponds

	Python call: Pond(xc,yc,rw,gamma,parent)

	Attributes:

	  - xc,yc: center of Pond (inherited from Circle)
	  - radius: radius of Pond (inherited from Circle)
	  - parent: parent of the element (inherited from Element)
	  - parameters: ColumnList of strength parameters, set to gamma, positive for exfiltration 
	  - solver: solver (inherited from Element)
	  - elementList: list of elements, defaults to itself (inherited from Element)
	  - __numberOfParameters: number of parameters, set to 1

	Overloaded functions from Element: __init__, potentialInfluence, dischargeInfluence
	'''

	def __init__(self,xc=0.0,yc=0.0,rw=0.0,gamma=0.0,parent=None):
		'''Constructor'''
		# Call base class constructors
		geometry.two_d.Circle.__init__(self,xc,yc,rw)
		base.Element.__init__(self,parent)
		self.__numberOfParameters = 1
		self.parameters = ColumnList([ gamma ])

	def __repr__(self):
		# Parent not yet included
		return 'Pond' + str((self.center.x,self.center.y,self.radius,self.parameters[0]))	

	def __complexPotentialInfluence(self,xy):
		if self.isInArea(xy):
			return complex(0.25 * ( xy.distance2(self.center) -
					        self.radius*self.radius ), 0)
		else:
			return 0.5 * self.radius * self.radius * \
			       log((complex(xy) - self.zc)/self.radius)
	def potentialInfluence(self,xy,t=0):
		assert isinstance(xy,XY),"XY point is required"
		return RowList([(self.__complexPotentialInfluence(xy)).real])

	def __complexDischargeInfluence(self,xy):
		if self.isInArea(xy):
			return -0.5 * complex(xy.x-self.center.x, -(xy.y-self.center.y))
		else:
			return -0.5 * self.radius*self.radius / (complex(xy) - self.zc)

	def dischargeInfluence(self,xy,t=0):
		assert isinstance(xy,XY),"XY point is required"
		comdis = self.__complexDischargeInfluence(xy)
		return RowList([XY(comdis.real,-comdis.imag)])

	def rechargeInfluence(self,xy,t=0):
		assert isinstance(xy,XY),"XY point is required"
		if self.isInArea(xy):
			return RowList([-1.0])  # Since recharge function returns what goes in
		else:
			return RowList([0])

	def nearElement(self,xyz1,xyz2,step):
		xy1 = xyz1.XY()
		xy2 = xyz2.XY()
		rv = TraceCheck()
		d1sq = self.distance2(xy1); d1 = sqrt(d1sq)
		d2sq = self.distance2(xy2); d2 = sqrt(d2sq)
		ssq = xy1.distance2(xy2); s = sqrt(ssq)
		oneplus = 1.000001
		# If point 1 outside pond and point 2 inside pond, or vice versa, step onto edge
		if d1sq > self.radsq and d2sq <= self.radsq:
			rv.change = 1
			cost = (d1sq + ssq - d2sq) / (2 * d1 * s)
			snew = ( 2*d1*cost - sqrt(4*d1sq*cost**2 - 4*(d1sq - self.radius **2)) )/ 2
			rv.xyzchanged = xyz1 + oneplus*snew/s*(xyz2-xyz1)
		if d1sq <= self.radsq and d2sq > self.radsq:
			rv.change = 1
			cost = (d2sq + ssq - d1sq) / (2 * d2 * s)
			snew = ( 2*d2*cost - sqrt(4*d2sq*cost**2 - 4*(d2sq - self.radius **2)) )/ 2
			rv.xyzchanged = xyz1 + oneplus*(s-snew)/s*(xyz2-xyz1)
		return rv
